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1.0 Summary 


“This research and development effort has produced the following 


* Software suitable for the construction of finite element 


gravity disturbance fields. 


Software suitable for calculation of the gravity distur- 
bances within finite element fields: 
- using Chebyshev polynomials 


- using Orthonormal polynomials 


Example finite element gravity disturbance fields have 
been developed and tested. Relevant comparisons bet- 
ween the Chebyshev polynomials approach and the Ortho- 


normal polynomials approach have been made. 


Optimizations have been incorporated to increase compu- 
tational efficiency and to determine a model for the 


geopotential, using Orthogonal polynomials. 
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2.0 Preface 4 
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This report constitutes the final report under Contract No DAAg70- 


. 78-c-0072 pérformed by the Virginia Polytechnic Institute and State 
~Wtiveretty for the U. S. Army Engineer Topographic Laboratories, Fort 
Belvoir, Virginia, under the sponsorship of the Defense Mapping Agency. 

' ‘The current effort is primarily concerned with the investigation of 
new choices of computationally more efficient basis functions which can 
be used to model the local fine structure of the gravity vector. This 
development, to date, has produced two separate versions of software 
systems. This report documents each of these systems as well as their 
theoretical background. Comparisons are made based on the result of 
several test cases. Also included are detailed discussions of the 
software and standard print-outs of each version. In addition, the 
programs contain enough detailed comment statements to allow the know- 
ledgeable user to use and modify particular sections to suit his or her 
own needs. Finally, the report also indicates a possible approach to 
determine a model for the geopotential itself, which at the same time 
would increase the computational efficiency of determining the gravity 
vector. 

The authors very much appreciate the technical liaison of Mr. L. A. 
Gambino, Director, Computer Sciences Laboratory (USAETL), who served as 
technical monitor of this work. Discussions with Dr. R. W. Ballew, Dr. 
B. Louis Decker, and Mr. H. W. Howard (all of the Defense Mapping Agency, 
Aerospace Center) were also helpful in defining the objectives of our 
effort. The excellent programming support of Mr. John T. Saunders and 


Mr. John J. Smith is gratefully acknowledged. 
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3.0 Introduction 


The central theme in the development of the present software sys- 
tems is the finite element approach to piecewise approximation. The 
finite element concept may be considered as an extension to three dimen- 
sions of one dimensional piecewise approximation techniques. Qualita- 
tively, some quantity that is functionally related in a "complicated" 
fashion to position over a "large" interval may be replaced for much 
"smaller" intervals of position by much "less complicated" approximating 
functions. Thus the complicated function, which will normally be "more 
expensive" to evaluate, will have been replaced by a set of locally 
valid functions which will normally be "less expensive" to evaluate. 

The idea may easily be extended to three dimensions. A quantity 
which is functionally related to position over some large region may be 
replaced by a set of approximations, each of which is valid only ina 
small local volume; the set of local volumes spanning the large region 
of interest. The local approximations are typically much easier (faster 
computationally) to evaluate. The finite element concept may be applied 
to any of a variety of modeling problems; this investigation centered 
specifically on fine structure gravity modeling. 

Previously, earth-fixed spherical coordinate based finite element 
fields have been successfully used to replace globally valid spherical 
harmonic representations of the geopotential and its derivatives. But 
because the model (in the present study) used to simulate gravity 


disturbance data (Model 310) consists of a set of point masses placed in 


relation to the Geodetic Reference Surface of 1967, ellipsoidal coordinates 


(H, A, ¢ - see Figure 1) were used in the finite element modeling process. 


Finite element fields may be defined by ellipsoidal coordinates as 


shown in Figure 2. The region to be modeled is defined by a set of 


maximum bounds (H * 
ma 


> dnax? We) and a set of minimum bounds (Hoan? r 


x min 


). The region is then divided into a set of smaller volumes, or 


min 
finite element cells, where each cell has dimensions H 5 x ; 
cell’ cell 
rt) . Thus the 300 km x 10° x 10° region above the reference ellipsoid, 
cell 


altitude 0 - 300 kn, 
longitude 70°E - 80°E, and 
latitude 25°s - 35°s, 
could be modeled by one hundred (100) cells each 300 km x ae ae or 
) 


four (4) cells each 300 km x 5° x 5° or one cell 300 km x 10° x 10°, 


etc. 


The approximation valid in any one cell is generated independently 
of the approximation for any other cell. A set of gravity disturbance 
“observations” are produced for each cell directly from Model 310 and 
then approximated for only that cell in a least squares fitting process. 
The approximating models thus generated are then stored on a high speed 
rotating mass storage device for later use in the computation of the 
approximate gravity disturbance in any finite element cell of the region 


modeled. 
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Figure 1 Earth-Fixed Rectangular and Ellipsoidal 
Coordinate Systems. 
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Figure 2 Ellipsoidal Finite Element Field. 


os 4.0 Finite Element Formulation 


ij Define the gravity disturbance, dg: 
is Sactual . reference 6g (1) 
Cs } 
| th lon + Sactual ~ Breference (2) 
eee 


Now, replace the components of 6g with locally valid polynomial approxi- 


mations, 6G, of the form, 


va “uP sk 
N n n-i 
6G = wane = soy ate “tal £ je %1 »X2 0X3) 
SG orth “north, 5 
(3) 
where, 
N = NORDER = the highest order polynomial, 
k=n-i- j 
*% eo Bee! iced 
X, = (QQ - Nef cell (non-dimensional coordinates) (4) 


tal 
i] 


3 = ret /* cell 


fiik = chosen polynomial basic functions, and 


Cc » C » and C are constant coefficients determined via 
UPS 5k sais aie | 
™m m 
ie 2 ey 2 
least square fits so that a (66, 88up) 4 ; 2 a SE | Fe Paar 


m 
and jy (86 orth - S8north i2 are minimized for some local volume of 


space (i.e. the cell being modeled). 


Thus, to determine the three sets of constant coefficients, Cup? 


Cc » and C for each cell, the following linear least squares 


east north’ 


problems must be solved, 


[166 = SBup/ | = [AC - Sg y II = min. 
[I8Scase ~ SBease!! ~ |lAScage ~ Beast! = min. ) 
[166 orth = SBorth | Fr [Ac orth ¢ SSoorth! | +e 
where the coefficient matrix, A, has the form, 
Foo0*? — fgg1? Foro) ds Sung 


Foo0'2) = F991 (2) fgg (2) s+ Egg (1) 
fm (6) 
Foo0™ — F991 — Foyg(M +++ Engg 


in which m is the number of observations in the cell and fii 


represents f,,,(X; ,X2, X3) corresponding to the th observation. 
oe 


Furthermore, the three coefficient vectors, £55? Cae? and Seth? 
have the form : 
C C Cc 
UPogg east yon northoog 
c C 
UPog1 east) “northyy, ’ 
Cc 
“Po10 “eastoyo “northy 9 
Lip ={ “Pioo )» Seast 88%100 )» Snorth Bere Oy 
Cc Cc 
UPog2 easton “northy, 
Cc Cc C 
“Poi eastyi) north), 
Cc C Cc r 
UPN eastyog northyog 4 


where i 
N = NORDER 


and the three sets of gravity disturbance observations, a | 


and Of eh? are 


Sek, ok, 45, ) ég Se Sy Se be, 
up 1, 24 34 east 1, 21 35 
68,0 (x Ri ae $Beast 1, 1%2,°%3,) 
6 = 6 
Eup sa aes uae Beast = so ae Ms 
a es Se *Sease 1 72 *%3) 


bg (> Sr? ae oa | 
north 1, 25 35 (8) 
Se orth 2X, »X, ) 
2 2 2 
6 _— 
north = Saaeeh 14/8213.) 
SB porth 1 »Xo ie ) 
m “m ~m 
where a i=1, 2, ..., m} = a specified set of (H;» die ¢,) 


coordinates - usually a uniform observation grid in any one cell. The 
gravity disturbance observations are determined directly from Model 310 
evaluations on the observation grid in each cell. 

Following the standard least squares approximation procedure we ob- 


tain, for the “up-direction" 


Cc =B "Sup (9) 


~up 


B= (a! a) ta? 


Similar results are valid in the “east- and north-directions". 


5.0 Chebyshev Polynomials Finite Element Models 


The choice of basis functions in Eq. (3) is still open. Several 
sets of basis functions have been used to date. Polynomial sets are of 
particular interest because of the ease with which they can be manipulated. 
The most obvious set is given by 


= Yeviye by 
Es jx Bp -%poXq) = X)X5X, i,j,k = 0,1,...,N (11) 


Another set consists of products of shifted Chebyshev polynomials 


fs je Ky oXqoXqo) = TYCK)T, (HTK) 5, = 0,1,..65N (12) 


j 


Chebyshev polynomials, t (x); of order n may be generated as follows: 


PON SR SNR ENTS NERA OT INN 


n t (x) 
0 1 
a : | x 

2 2x2 - 1 


For n > 2, a recursion formula may be used: 
tO) * 2X ty = t 26%)» =i Xesek (13) 


Shifted Chebyshev polynomials, T(x) may be computed by substituting 


X= 2X-1 
so that, 
TT, = t (x) = t, (2x - 1), Q< kX <1 (14) 


To generate shifted Chebyshev polynomials a set of X's with values bet- 
ween 0 and 1 are needed. The difference between some point (H,A,$) in 


il} w 
a cell and a reference point (CH eve? reg? (the "lowermost corner 


of the cell) is divided by the cell's dimensions, r 


Aeei1’ *ce11’ *%ce11’ 
as in (4) to obtain a set of non-dimensional coordinates, X)> Xo» X,, 
for each cell. The non-dimensional coordinates may then be used to gen- 


erate the shifted Chebyshev polynomials. 


5.1 Discussion 


Certain aspects of the process by which finite element fields are 
created may be capitalized on under special circumstances. When it is 
known beforehand that uniformly gridded gravity disturbance observations 
will be available, only one least squares matrix (A matrix) need be 
generated. This is true because the set of non-dimensional coordinates, 


{(X, »Xy »K, ), i=1, 2, ..., m} of one cell's observation grid will be 
i i 5 i 


the same as every other cell's non-dimensional observation grid coordinates 
when the observation grid pattern is set for all cells. If the same 
number of observations are made in each cell and the positions of the 
observations points are the same relative to the reference point (H 


ref’ 


r ) in each cell, then all cells will have the same A matrix. 


ref’ * ref 
The A matrix, after all, consists only of products of the locally evaluated 
polynomials, which are functions only of the non-dimensional coordinates 
of the observation points. 

This all means that only one A-matrix need be generated for the 
first cell of a finite element field; but the gravity observations must 
be generated for each cell. A reduction of the A-matrix to upper tri- 
angular form may be performed for the least squares fitting process, 


after which each set of observations may be similarly reduced and then 


back-substituted to produce a set of coefficients for each cell. The A- 


matrix, again, needs to be reduced only once. 


It should be noted that this method of finite element field generation 
does not rely on Model 310 or any other particular gravity disturbance 


model. Model 310 could be replaced by any other process that is capable 


os of producing the gravity disturbance observations at the observation 


grid points of each cell. 
Once a finite element field has been created, the Chebyshev poly- 
nomial gravity disturbance approximation 6G may be computed by evaluating 


(3). The subscripts ijk always follow the same pattern for a given 


is NORDER. This means that a “current coefficient number", 2, may be 


attached to each group of subscripts ijk and that the total number of 
coefficients, NC, necessary to model any one component of the gravity 
disturbance in one cell will be constant for a given NORDER. Table 1 
shows these relationships for various NORDER's. 

Since we can precompute and save the subscripts, ijk, in the appro- 
priate pattern such that they are functions of only the current coefficient 


number, £, (3) can be written, 


8G, Cs 
NC - 
a a 3 Coast, (M4 6Xq)T, Kp) T, Ky) (15) 
k=1 < 
SC orth & 


northp 


where, 


NC = the total number of coefficients for one component of the 
gravity disturbance, as above 


ijk = precomputed and saved as functions of 2 


Xj »X_sX3 = functions of H,A,$ and the constants Heo? def? ef? 
Heel? dell? el] as in (4) above. The total number of summations in 
(15) is the same as in (3) but the amount of bookkeeping has been cut 
down tremendously by doing business in this fashion. Computational 
savings of up to 50% were obtained by using (15) in lieu of (3) for 
evaluation purposes. Table 1 was also used in the generation of the A- 


Matrix, but the saving here is not as significant due to the previously 


discussed fact that the A-matrix is only generated once. 
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NORDER ijk if NC 
1 


0 000 1 
1 ool 2 
010 3 
100 4 4 
2 002 5 
oll 6 
020 7 
101 8 
110 9 
200 10 10 
3 003 11 
012 12 
300 20 20 
N 0 N 
01 N-1 
: 3 M1 oon 
N00 NC } Jao 
n=l m=1 


Table 1 Subscript Patterns for Various NORDER's. 


6.0. Orthonormal Polynomials Finite Element Models 


As demonstrated by the numerical results, Chebyshev polynomials 
yield an entirely satisfactory solution. However, by choosing a judicious 
set of orthonormal polynomials as basis functions, it is possible to 
avoid the computation of the matrix inverse in Eq. (10). In addition, 


it will be possible to compute the coefficients c independently of 


ijk 
each other. This is especially convenient when the order of approximation 
is increased for no need exists to recompute the least squares operator 


B in Eq. (10), or any lower order coefficients. 


Let us introduce the following basis functions 
£5 51% 2% 9X5) = £ (Kf; (XE, (X3) i,j,k = 0,1,...,N (16) 


Where f,, f fi. represent (to be specified) one-dimensional polynomials 


i’? Sag 
in X)> Xo» and X, respectively. 
Using Eq. (6), it is not difficult to show that the matrix Ala has 


elements of the form 


m 
Y Specht Oe). gage Oylyee eh 
ia IE SRY (17) 
a sB,y = 0,1,...,N 
and introducing Eq. (16) into Eq. (17) yields 
~ my my 
r L £,(2,)£,(2,) P Lo f, (Lyf (Lp) | Lo £.(L3)£ (25) (18) 
1 2 3 


where m, (i = 1,2,3) represents the number of observations in the X,- 
direction, and m is the total number of observations in the cell. 
Next, let us assume that the set £, (X)) can be constructed ortho- 


normal over the discrete domain of interest. Hence, 


») 
Fete yee) © 8 (19) 
£,=1 
1 
Similar expressions can be written for the sets £5 (X,) and £, (X3)- 


Substituting Eq. (19) into Eq. (18) yields 


85455 e°Ky (20) 
and because (20) is unity if and only if i=a, j = 8, k = y, we conclude 
that the matrix Ata reduces to the unit matrix. 
The problem is now reduced to the construction of an orthonormal 
set of one-dimensional polynomials £, (x) (i = 0,1,...,N). 


Any orthogonal set of polynomials has the following three-term re- 


currence relation[3] 


i 6c Ait i) * arurest eta) = Bee ee) (21) 
Ry, itl i @, yh, ‘i-1 4 


where hy stands for the leading coefficient in £ , (x) and By is an unknown 


constant, and where we use the inner product notation 
2 <f,,8> = j £7(x)dx (22) 
i 5 idea 8 Di 


and where D is the domain of interest. If we choose the following non- 
dimensional coordinates 


y 2(H-H_) 
Ee Beas 


20-2) 
Xy aie a (23) 
cell 


2(6-6,) 
%ee1l 


3 


then, the domain D is (-1,1) which enables us to conclude that the 


functions fo = 1 and fy = x are orthogonal over D. Moreover, using the 


Gram-Schmidt procedure it can be shown that fo, (k = 0,1,...) is even 


while fot (k = 0,1,...) is odd. From Eq. (21) it then follows that all 


B, must be equal to zero. Finally, because ho = hy = 1, we also must 


have hy = 1 (i = 2,3,...). Relation (21) then reduces to 


a ig 


<f;_,f 


f (x) = xf, - 


i+] a (4) 


a" 


which, together with fo = 1 and fi = x, yields a recurrence relation for 


a set of orthogonal polynomials. Normalization is accomplished by 


i. wees ie Bb cea (25) 
i ,,£ 17? 


a 


where fy represents the normalized polynomial. For a discrete observation 


grid we replace Eq. (22) by 
Y #0) (26) 


It is also necessary that the local coordinates of the grid points in a 
particular direction add up to zero. Furthermore, we must have that the 
number of grid points in a particular direction is at least Ml. | 

Specifically Table 2 gives the orthonormal polynomials £,G = 0,1,2,3) 


for the grid points (-1,-1/3,1/3,1). } 


18 


fy = +500000 
f, = +670820x 

£, = 1.125000x" - .625000 
f, = 2.51558X° - 2.29197X 


Table 2 Orthonormal Polynomials 


Introducing the polynomials (25) into Eq. (16) and into Eq. (6) 
reduces Eq. (10) to 


es (27) 


so that, from Eq. (9) it follows that 


ae 
Se Bey (28) 


It is then possible to write an explicit expression for each coefficient 


as follows 


ijk 
m 
“15k, ho Es 5460) 58,,, (2) (29) 


Note that each coefficient c can be computed independently and no 


as as 
need exists to recompute B when the value of N changes. This feature is 
very valuable when, for example, a prototype gravity model of a launch- 
region is desired but it is difficult to specify the order N required 
for an accurate representation; also note that N could be variable de- 


pending upon the frequency and amplitude of local gravity anomalies. 


Also, in case of a global gravity model, we can use a truncation of a 


spherical harmonic series for 82, and compute 


“ik, = Ip £5 5x°8up (30) 


by carring out analytical integrations. This reduces Eq. (3) to a 
truncated Fourier series. Because a continuous spectrum of "measurements" 


is used, we can expect more accurate approximations. 


6.1 Discussion 


The recursive relationship of the orthonormal polynomials is very 
similar to the relationship previously used with Chebyshev polynomials. 
This formulation however does contain two additional division operations 
as well as a square root evaluation which must be performed on each 
polynomial. If one division and the square root values are saved during 
the preliminary calculation (SUBROUTINE ORTHO) and passed on for the 
actual evaluation, (SUBROUTINE MULT) a time savings can be realized at 
the small expense of a few (42) additional "words" of storage. 

Note that even with this saving, the original Chebyshev polynomials 
can still be evaluated slightly faster, the difference only aitiss the 
execution of a single division operation. However, when the advantages 
noted earlier are taken into account, the choice of orthonormal poly- 
nomials may be more beneficial, especially considering the effort that 
is involved with several different "set-ups" of variable polynomial 
orders and/or observation grid patterns where the Chebyshev version 
would require extensive computational time to process the least squares 
matrix and evaluate the coefficients. Indeed, the orthonormal approach 
does not require evaluation of (atay~t and in addition, every coefficient 


os 5k is independently computed. 


a Al i ah lh 


7.0 Discussion of the Software 


A variety of Fortran programs and subroutines were written during 
the course of this research effort. The software developed by Junkins 
and Remdeta’ was converted from its original CDC form to run on the 
VPI&SU IBM System 360 Model with 158 Processors. After correctly icpro- 
ducing their test results, the software was modified to incorporate some 
programming efficiencies. From this starting point, new programs were 
developed to implement the orthonormal approach as well as test various 
additional approaches. 

In general, the final software package can be divided into three 


separate stand-alone modules, regardless of the basic functions chosen. 


Section I: A Mass Model 310 Initialization Section. 
q Section II: A Least Squares Reduction and Coefficient Determination Section. 


Section III: A Test Section. 


Section I, A Mass Model 310 Initialization Section, consists of a 
id single fortran program, MASPOS, and its coded input data. This job step 
TT is identical for all the versions of the software developed so far and 
it will be required by any new version which makes use of Mass Model 310 
to supply simulated gravity disturbances for use as observations or as a 
control to compare the modeling equation's calculations against. Clearly, 
this module can be replaced by input software if other sources (e.g., 
measurements) of gravity disturbances are available. 


This program presently inputs coded values of the 1080 point masses 
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of Mass Model 310. Their location's are converted from the geodetic 


coordinates (H,A,¢) to rectangular coordinates (x,y,z) and the mass 
codes are converted to their actual values. This information is then 
written onto a sequential access file (FILE 2), for use by SUBROUTINE 
PTMASS, which is described elsewhere. More specifically, PROGRAM MASPOS 
creates an unformatted sequential file, FILE 2, which contains the 
precomputed products of the gravitational constant and each point mass 
of Model 310 and the precomputed geocentric rectangular coordinates of 
each point mass. PROGRAM MASPOS requires as input coded information 
about Model 310 in the form of 30 cards. Each card represents 36 mass 
points in a grid row of equal latitude where the mass points are 50 
longitude minutes apart. The mass points should be listed on each card 
by increasing longitude and the cards should be arranged in order of 
decreasing latitude. PROGRAM MASPOS reads only the coded multiplication 
factors for each of the 36 mass points of the grid row, where the coded 
factors indicate the mass of the respective mass points, 1 = -1.x10!9 
grams, 2 = 0. grams, and 3 = +1.x10!9 grams. The coded factors are 
located in columns 25-60 of each data card. All other data card columns 
are ignored (for information only, columns 1-9 indicate south latitude 
in arc minutes of the grid row and columns 10-18 indicate east longitude 
in arc minutes of the west-most mass point in the row). 

Section II, A Least Squares Reduction and Coefficient Determination 
Section, is where the actual modeling equation coefficients are calculated. 
Since the coefficients are determined by the choice of the basis functions 
and the choice of basis functions is a primary issue under study, each 
version of the software for Section II differs in several important aspects 


while the fundamental sections of the programs, such as establishment of 


the finite element cells and grid pattern and determination of the 
gravity disturbances will be the same. | 

The main program for Section IIA, LOCALG, performs all the neces- 
sary preliminary calculations that are required by each version. Using 
card input, the finite element field is established and broken down into 
finite element cells and the observation grid pattern is also esta- 
blished. The critical field information is then written onto a storage 
file (FILE 3) for later use by PROGRAM FINEG. More specifically, 
PROGRAM LOCALG generates the A-matrix for a "typical" finite element 
cell and a set of gravity disturbance observations for each cell of a 
finite element gravity disturbance field. The gravity disturbance 
observations are evaluations of SUBROUTINE PTMASS on the uniform ob- 
servation grid in each cell. The user defines a finite element gravity 
disturbance field by providing the following items as inputs to PROGRAM 


LOCALG: 


1. The integral order of the locally valid polynomial 
approximating functions, NORDER, (0 < NORDER < 6). 

2. The observation grid pattern (same for all cells) 
to include the number of observations in the up, 
eastern and northern directions, MOBSU, MOBSE, 
MOBSN. The total number of observations (MOBSU x 
MOBSE x MOBSN) should be at least three times the 


number of coefficients, NC (see Table 1), in the 


approximating model, Furthermore, the number of 


These thirteen items are input on four data cards. 


observations in each direction should be one 

more than NORDER (MOBSU = MOBSE = MOBSN = 

NORDER + 1). The maximum number of observa- 

tions allowed by the present program dimensions 

is 343 (a 7x7x7 observation grid). Non-uniform 
sample grids (e.g. 5x7x7, 6x5x4, etc.) have been 
found to produce unpredictable results - often 

an unacceptable fit. 

The size of each finite element cell in H,A,$¢ 

must be specified, HCELL, ALCELL, APCELL. HCELL 

is cell size in H (meters), (HCELL < [HMAX - HMIN]). 
ALCELL is cell size in A (degrees), (ALCELL < 

[ALMAX - ALMIN]). APCELL is cell size in » (degrees), 
(APCELL < [APMAX - APMIN]). 

The minimum and maximum H,i,> boundaries of the 
finite element field, HMIN, ALMIN, APMIN, HMAX, 
ALMAX, and APMAX. HMIN and HMAX are heights (H) 
above the reference ellipsoid in meters, (HMIN < 

O and HMAX > HMIN). ALMIN and ALMAX are ellipsoidal 
longitudes (A) in degrees, (-180° < ALMIN < ALMAX 


< 180°). APMIN and APMAX are ellipsoidal latitudes 


(9) in degrees, (-90° < APMIN < APMAX < 90°), 


MOBSE, and MOBSN are input in (4110) format on the first card. The 


remaining three cards must have a (3E20.14) format, with the finite 
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NORDER, MOBSU, 


element field delimiters, HCELL, ALCELL, and APCELL on the second 
card, HMIN, ALMIN, and APMIN on the third card, and HMAX, ALMAX, and 
APMAX on the last card. 

Because gravity disturbance observations are generated within 
PROGRAM LOCALG by SUBROUTINE PTMASS, PROGRAM LOCALG also requires as 
input the precomputed products of the gravitational constant and each of 
the mass points of Model 310 and the precomputed rectangular coordinates 
of each mass point. These quantities should have previously been stored 
sequentially via an unformatted write on FILE 2. PROGRAM LOCALG accesses 
these quantities with an unformatted read of FILE 2. 

PROGRAM LOCALG produces as output an unformatted sequential file, 
FILE 3, which contains the A-matrix for a typical finite element cell 
and a set of gravity disturbance observations for each cell of the user 
specified finite element field. This file is properly formatted for 
later use by PROGRAM FINEG. PROGRAM LOCALG requires two supporting 
subroutines - SUBROUTINE CHEBY (or SUBROUTINE ORTHO) and SUBROUTINE 
PTMASS. 

The basis polynomials are evaluated by the subroutines CHEBY or 
ORTHO depending on the particular version, and their outputs are used to 
construct the least squares matrix. 

In particular, SUBROUTINE CHEBY (X, N, TA) returns Chebyshev 
polynomials in array TA, evaluated at a specified value, X, through the 
specified order, N. TA(1l) = Ty (X), TA(2) = T,@), TA(3) = T, (x), etc. 
TA is a seven element array; hence the specified order is bounded, 

0 < N < 6. SUBROUTINE ORTHO together with SUBROUTINE MULT perform the 


same function in the case of orthonormal polynomials. 


In all versions, the gravity disturbance observations are det- 
ermined for each cell by SUBROUTINE PTMASS, the final subroutine in 
Section IIA. 

SUBROUTINE PTMASS is used both to produce observations for the 
least squares fitting process and to produce observations to test 
against for error analyses of finite element fields. Because SUBROUTINE 
PTMASS required the precomputed products of each point mass and the 
gravitational constant and the earth-fixed Cartesian coordinates of each 
point mass of Model 310, another small program, PROGRAM MASPOS, was 
necessary to calculate and store these quantities. SUBROUTINE PTMASS 
uses Mass Model 310 to determine the gravity disturbance observations at 
each cell grid point. 

The rectangular coordinates (x,y,z) of the point are sent to the 


subroutine and the gravity disturbances are determined by: 


5g 
x 
1080 
bg = Me y  PMVALS Ss (30) 
y =1 d, 
~ 6g ae 


where PMVALS ; is the actual value of point mass i which was calculated 


by MASPOS in Section I 


*,) x - POSITS, 
1 
di =< dy, >= <y - POSITS (31) 
o i, 
dz; zZ- igs 


POSITS , are the (x,y,z) coordinates of point mass i which was cal- 
n 
culated by MASPOS in Section I and 
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PROGRAM LOCALG LOGIC FLOW (CHEBYSHEV VERSION) 


START 
INITIALIZE CONSTANTS 


POINT MASS VALUES AND 
<= 


PRECOMPUTED POSITIONS FILE 2 

INPUT POLYNOMIAL ORDER, NORDER, MOBSU, MOBSE, 
OBSERVATION GRID PATTERN, MOBSN, HCELL, ALCELL, 
FINITE ELEMENT CELL SIZE, APCELL, HMIN, ALMIN, 

AND FINITE ELEMENT FIELD <== 6 APMIN, HMAX, ALMAX, APMAX 


BOUNDS 


GENERATE ALL THE BASIS 
FUNCTION VALUES AT EACH 
OBSERVATION GRID POINT 


ITE CRITICAL FINITE 
ELEMENT FIELD INFORMATION 
ONTO THE SEQUENTIAL ACCESS 
FILE 


===> FILE 3 


GENERATE THE LEAST SQUARES 
MATRIX AND WRITE IT ONTO ===> FILE 3 
THE SEQUENTIAL ACCESS FILE 


LOOP THROUGH THE FINITE 
ELEMENT FIELD, ONE CELL 
PER LOOP 


GENERATE THE GRAVITY 
ANOMALY OBSERVATIONS 


(S815 SBeast 08 
WLTH PTMASS 


north) 


WRITE THE THREE OBSER- = FILE 3 


VATION VECTORS ONTO THE 
SEQUENTIAL ACCESS FILE 


GO TO THE NEXT FINITE 
ELEMENT CELL BY FIRST 
VARYING $ THEN i AND 
FINALLY H UNTIL ALL 
THE CELLS HAVE BEEN 
OBSERVED 


STOP 


PROGRAM LOCALG LOGIC FLOW (ORTHONORMAL VERSION) | 


START 
INITIALIZE CONSTANTS 


READ POINT MASS VALUES 


—=—=—=—=————= 3 ; 
AND PRECOMPUTED POSITIONS < FILE 2 


INPUT POLYNOMIAL ORDER, NORDER, MOBSU, MOBSE, 
OBSERVATION GRID PATTERN, MOBSN, HCELL, ALCELL, 
FINITE ELEMENT CELL SIZE APCELL, HMIN, ALMIN, 
AND FINITE ELEMENT FIELD SS. APMIN, HMAX, ALMAX, 
BOUNDS APMAX 


GENERATE THE COEFFICIENTS 
FOR THE BASIS FUNCTION IN 
THE H DIRECTION 


GENERATE ALL THE BASIS 
FUNCTION VALUES FOR EACH 
CONSTANT H PLANE 


GENERATE THE COEFFICIENTS 
j FOR THE BASIC FUNCTION IN 
a THE X DIRECTION 


GENERATE ALL THE BASIS 
FUNCTION VALUES FOR EACH 
CONSTANT A PLANE 


% GENERATE THE COEFFICIENTS 
4 FOR THE BASIS FUNCTION IN 
i THE. ¢ DIRECTION 


GENERATE ALL THE BASIS 
FUNCTION VALUES FOR EACH 
CONSTANT $ PLANE | 


WRITE THE CRITICAL FINITE 
ELEMENT FIELD INFORMATION AND FILE 3 
ORTHONORMAL COEFFICIENTS ONTO => 


THE_SEQUENTIAL ACCESS FILE 


GENERATE THE LEAST SQUARES MATRIX 
AND WRITE IT ONTO THE SEQUENTIAL —= FILE 3 
ACCESS FILE 


v 


FIGURE 4 
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PROGRAM LOCALG LOGIC FLOW (ORTHONORMAL VERSION) CON'T 


LOOP THROUGH THE FINITE 
ELEMENT FIELD, ONE CELL 
PER LOOP 


GENERATE THE GRAVITY 
ANOMALY OBSERVATIONS 
(58,4? Beast Snorth? 

WITH PTMASS 

WRITE THE THREE OBSERVATION 

VECTORS ONTO THE SEQUENTIAL == FILE 3 
ACCESS FILE 


GO TO THE NEXT FINITE ELE- 
MENT CELL BY FIRST VARYING 

¢ THEN A AND FINALLY H UNTIL 
ALL THE CELLS HAVE BEEN 
OBSERVED 


STOP 


FIGURE 4 CON'T 


2 2 2 
dq, = dx, + dy; + dz, 


In particular, SUBROUTINE PTMASS produces gravity disturbance "obser- 


vations", 88,59 6g 


mane SB orth? to be used by PROGRAM LOCALG and 
PROGRAM FINEG for fitting purposes and to be used to test as data 
against SUBROUTINE FINITE. SUBROUTINE PTMASS produces the gravity 
disturbance data by using the 1080 point masses of Model 310. SUB- 
ROUTINE PTMASS requires the COMMON/MASPOS/ be filled with the precom- 
puted products of the gravitational constant and each point mass and the 
geocentric rectangular coordinates of each point mass. SUBROUTINE 
PTMASS interfaces with the calling program via COMMON/XYZ/ through which 
the rectangular coordinates, x, y, z, of the calling point are received 
and the components of the gravity disturbance og.» = and 58, (DELGX, 


DELGY, DELGZ) are transmitted. The rectangular components of the gravity 


t? and 


disturbance must be transformed to geodetic components BE oc? OB aa 


$8 orth’ outside SUBROUTINE PTMASS (i.e., in the calling routine). The 
logic flow of PROGRAM LOCALG for the different versions are illustrated 
in Figures 3 and 4. 

The Chebyshev and Orthonormal versions both use a second section, 
Section IIB to actually calculate the modeling coefficients. These two 
versions are different enough to warrant separate discussions. 

In the Chebyshev version, the critical field information and the 
least squares matrix are read in. This version must calculate the 
modeling coefficients by solving the least squares solution. In lieu of 
solving the classical normal equations, a more efficient matrix re- 


duction algorithm (ALSQ) is employed. A sequence of Householder ele- 


mentary matrices are determined which, when premultiplied into the A-matrix 


= 


reduce it to a very particular upper triangular form. This triangular 


system can then be solved by backward substitution for the coefficients. 

The A-matrix and its dimensions are passed to the main entry point 
of SUBROUTINE ALSQ. This section of ALSQ reduces the least squares 
matrix to its upper triangular form which is returned to FINEG. Next 
the gravity disturbance component vectors of each cell are read in one 
at a time. The reduced least squares matrix and the observation vector 
are then passed to the ALSQ1 entry point uf SUBROUTINE ALSQ. This sec- 
tion operates on the observation vectors with the same sequence of 
Householder transformations, then solves the triangular system for the 
vectors of coefficients. The process is repeated for each of the up, 
east, and north measurement data to determine the corresponding coeffi- 
cient vectors. The coefficients of each cell are written into separate 
records of a random-access file (FILE 3) to be used by SUBROUTINE FINITE 
in Section III. In particular, SUBROUTINE ALSQ (A,Y,B,R2,NN,MM,NA) 
solves the linear least squares problem, ||AB - Y|| = minimum. SUB- 
ROUTINE ALSQ requires as input the coefficient matrix, A, the observa- 
tion vector which is to fit, Y, the number of rows used in the A-matrix, 
NN, the number of columns used in the A-matrix, MM, and the first 
dimension of the A-matrix, NA. NA must be equal to one plus the maximum 
number of rows in the A-matrix (NA = NN ax + 1). SUBROUTINE ALSQ returns 
B, the coefficients of the fit and R2, the sum of the squares of the 
residuals. 

The Orthonormal version of Section IIB produces exactly the same 
result but the process is much simplier. The simplification is based on 


the results found in section 5.0 of this report which shows that the 
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least squares coefficients are given by Eq. (29). Thus, the only pro- 


cessing needed, is to multiply the transpose of the unmreduced least 
squares matrix (i.e., the A-matrix) by the observation vectors (one at a 
time) to produce the corresponding modeling coefficient vectors (one at 
a time). Since this was basically the job performed by the secondary 
entry point, ALSQ]1 in the Chebyshev version, the new subroutine for the 
Orthonormal version is called ALSQl. (It does not require a secondary 
entry point.) The coefficients are again returned to the main PROGRAM 
FINEG and written (in a set of three, one for each component) onto the 
appropriate record of the random access file (FILE 1). The logic flow 
charts of both of these versions are given in Figures 5 and 6. In 
particular, PROGRAM FINEG generates sets of locally valid coefficients 


ged » and C for each cell of a finite element field which was 


Cup east north 
previously specified by inputs to PROGRAM LOCALG. The only input required 
by PROGRAM FINEG is the unformatted sequential file, FILE 3, previously 
produced by PROGRAM LOCALG. PROGRAM FINEG produces as output a random 
access file, FILE 1, which contains, in effect, a finite element gravity 
disturbance field. The first record of the random access file contains 
certain important information about the finite element field (NORDER, 
boundaries, etc.). In addition to the first record, the random access 
file consists of one record of coefficients for each cell of the finite 
element field. Thus the number of records in the random access file 

will be equal to the number of cells in the finite element field plus 1. 
As indicated, PROGRAM FINEG requires one supporting subroutine - SUB- 


ROUTINE ALSQ. 
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PROGRAM FINEG LOGIC FLOW (CHEBYSHEV VERSION) 
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. - 
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MODELING EQUATION ONTO THE 
RANDOM ACCESS FILE 

GO TO NEXT CELL 


STOP 
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i 
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PER LOOP 


seiihiah-anthidteent tient na aT ee ee 


READ GRAVITY ANOMALY OB- 
SERVATIONS FOR THIS CELL 
FROM THE SEQUENTIAL FILE 


COMPUTE POLYNOMIAL COEFFI- 
k CIENTS OF THE FINITE ELEMENT 
4 MODEL WITH SUBROUTINE ALSQ1 
WRITE COEFFICIENTS OF THE 
MODELING EQUATION ONTO THE 
RANDOM ACCESS FILE 
GO TO NEXT CELL 


STOP 


FIGURE 6 
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Section III, Test Section, is the portion of the software system 
which evaluates the accuracy and speed of the new modeling approaches 
verses the Mass Model 310 technique. The main program of this section, 
FINITES, is the same for both versions of the software which have been 
developed thus far. The only changes that are necessary are the DIMENSION 
statements which declare the type and number of the critical field data 
that must be read from record number 1 (one) of the random access file 
(FILE 1), the READ statement which performs this operation and the 
common statements that transfer it to FINITE. Everything else is the 
same and the procedure is as follows. 

After the critical field data are read in, the information produced 
in Section I for Mass Model 310 is also read in. This information is 
used to calculate the "actual" gravity disturbance measurements at each 
test point so that the modeling equation's values at the same points can 
be compared. Next, the intervals between test points (ISTEPH, ISTEPL, 
ISTEPP) are read in from cards along with the minimum (HMIN, ALMIN, 
APMIN) coordinates and maximum coordinates (HMAX, ALMAX, APMAX) of the 
test region. These coordinates are ellipsoicil. After the validity of 
the test pattern has been determined, the error analysis is initialized. 
The ellipsoidal coordinates (H, ALAM, APHI) of each test point are then 
passed to SUBROUTINE FINITE which returns the three gravity disturbance 
components which the modeling equation has calculated. These values are 
saved and the (x,y,z) coordinates of the same test point are determined 
and sent to SUBROUTINE PTMASS. This subroutine is identical to the 


subroutine PTMASS discussed earlier. It returns the three gravity 


disturbance components which Mass Model 310 has calculated. The 


disturbance gravity components from PTMASS are first converted to their 
ellipsoidal values and the components of the two techniques are compared 
on a point-by-point basis; the error analysis is printed out (the means, 
standard deviations, and absolute value of maximum errors are printed 
out along with the average time required to calculate each point). 

As noted previously, the actual evaluation of the gravitational 
disturbances is performed by subroutine PTMASS. For the Chebyshev and 
Orthonormal versions, although the basis functions are different, the 
procedure is exactly the same. First the ellipsoidal coordinates (H, 
ALAM, APHI) of the test point are received from the main program in the 
common block labelled HLP, and the point is compared against the modeled 
region minimum (HMIN, ALMIN, APMIN) and maximum (HMAX, ALMAX, APMAX) to 
see if the test point lies within the region modeled in Section II. If 
not, the gravity disturbance is evaluated by a call of SUBROUTINE PTMASS, 
a message indicating this fact is printed and the PTMASS values are 
returned to FINTES. If the point does lie within the modeled region, 
the critical field information is used to determine which finite element 
cell contains the test point. When the cell is determined, the three 
sets of coefficients (CU,CE,CN) are read from the random access file 
(FILE 1). Next the coordinates of the test point is normalized with 


respect to the finite element cell, and these values (X, »X 


X,) are then 


2°? 
passed to SUBROUTINE CHEBY for the Chebyshev version or SUBROUTINE MULT 


for the Orthonormal version and the basis functions are evaluated at 
that point. These subroutines are exactly the same as the ones used in 
Section IIA, to evaluate the basis functions to determine the least 
Squares matrix. The gravity disturbance components are then determined 


as the sum of the products of the coefficients multiplied by the 


evaluated basis function terms. These values are then returned to 
FINTES for error analysis. 
In particular, PROGRAM FINTES is a simple error analysis routine. 
It compares the gravity disturbance as approximated by SUBROUTINE 
FINITE with the gravity disturbance as evaluated by SUBROUTINE PTMASS 
(Model 310). Residuals are printed for each point of a user specified 
grid, along with the mean of the residuals, the root mean square of the 
residuals, and the maximum absolute value of the residuals, for each 
component of the gravity disturbance. The following items are required 
as input to specify the sample grid: 
1. The number of samples in each direction, ISTEPH, 
ISTEPL, and ISTEPP. (ISTEPH > 0, ISTEPL > 0, and 
ISTEPP > 0). The total number of sample points 
on the grid will be the product of these three 
numbers. 
2. The minimum sample grid coordinate values ("the 
lower corner"), HMIN (meters), ALMIN (degrees), 
and APMIN (degrees) and the maximum sample grid 
coordinate values ("the upper corner"), HMAX 


(meters), ALMAX (degrees), and APMAX (degrees), 


(0 < HMIN < HMAX, -180° < ALMIN < ALMAX < 180°, 

-90° < APMIN < APMAX < 90°). 
These nine items are input on three data cards. ISTEPH, ISTEPL, and 
ISTEPP are input on the first data card in a (315) format. HMIN, ALMIN, 
and APMIN are input on the second card in (F10.0, 2F10.2) format and 


HMAX, ALMAX, and APMAX are input on the third card in the same format. 
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PROGRAM FINTES requires three supporting subroutines - SUBROUTINE 
FINITE, SUBROUTINE CHEBY or ORTHO (see discussion of these routines 
above), and SUBROUTINE PTMASS (also discussed above). 

Because PROGRAM FINTES uses SUBROUTINE PTMASS, it also requires as 
input the unformatted sequential file, FILE 2, containing the precomputed 
products of the gravitational constant and each point mass of Model 310 
and the precomputed rectangular coordinates of each point mass. 

SUBROUTINE FINITE directly replaces Model 310 in applications 
involving the gravity disturbance 6g. Given the ellipsoidal coordinates, 
H, ALAM, and APHI (H,A,o), of a specified point, SUBROUTINE FINITE 


) 


of the gravity disturbance approximation, 6G. SUBROUTINE FINITE determines 


returns the ellipsoidal components, GU, GE, and GN ioe OE? ea 
which set of finite element coefficients to use at the specified point, 
reads them into core, and computes 6G. SUBROUTINE FINITE requires 
SUBROUTINE CHEBY or ORTHO (described previously) to produce basis poly- 
nomials at the specified point. Whenever the specified point does not 
lie within the bounds of the finite element field currently accessed, 
SUBROUTINE FINITE calls SUBROUTINE PTMASS, writes an error message, and 
sets the error flag ISITIN to zero. 

Prior to the first call to SUBROUTINE FINITE, the random access 
file, FILE 1, must be made available to the calling program, and two 
flags, IFLAG and ISITIN, must be set to zero and included in a COMMON in 
both the calling program and SUBROUTINE FINITE. This could be COMMON 
/HLP/ if necessary (COMMON /HLP/ is the vehicle by which SUBROUTINE 
FINITE receives the calling coordinates and returns the components of 
6G), although COMMON /IMARK/ was used. 

The following remarks are in order when (1)C FORTRAN is used. The 


requirements of CDC FORTRAN dictated that the CDC utility subroutines, 


ed 
wtebeee 


subroutines, OPENMS, WRITMS, READMS, and CLOSMS, be used to access the 
random access file. Furthermore, these CDC utility routines require 
that the user establish an index array (IMARK, in this case) for the 
random access file. If an index array is required by a particular 
version of FORTRAN, the index array, IMARK, should be included in a 
COMMON in both the calling program and SUBROUTINE FINITE. The logic 


flow for SUBROUTINE FINITE is displayed in Figure 7. 


8.0 Tradeoff Studies 


Several tradeoff studies were conducted on the VPI&SU IBM System 
360 Model with 158 Processors. The central processor execution times of 
SUBROUTINE FINITE were compared to SUBROUTINE PTMASS (Model 310). 
Although it is really a price that is paid only once and a priori in the 
laboratory, the central processor execution times necessary to create 
finite element fields were examined. Consideration was given to the 
tradeoff: total number of finite element coefficients versus NORDER 
versus cell size, for finite element models of a given region. Lastly, 
error analyses for a variety of NORDER's and cell sizes were conducted; 
specifically, the maximum absolute error versus NORDER versus cell size 
tradeoff was examined. 

Table 3 shows the execution time test results for three models: 
Chebyshev, Orthonormal, and Pointmass. SUBROUTINE PTMASS has a nearly 
constant run time. The run time for SUBROUTINE FINITE, however, depends 
on NORDER, the basis polynomial order, and on whether or not a random 
access call is needed. A random access call is needed anytime the set 
of coefficients currently in core is not the set for the cell in which 
the calling coordinates lie; that is, the previous calling coordinates 


were in a different finite element cell than the current calling coor- 
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SUBROUTINE FINITE LOGIC FLOW 


CALL FINITE (speciFy H,,% COORDINATES OF SOME POINT) 
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FINITE ELEMENT 
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~ 


Z COEFFICIENTS ; DISK 
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FILE 1 14 
COMPUTE FINITE ELEMENT CELL BOUNDARIES; NON-DIMENSIONALIZE 
THE CALLING COORDINATES; GENERATE BASIS POLYNOMIALS 


Compute 6G 
(MATRIX PRODUCT OF FINITE ELEMENT MODEL COEFFICIENTS 
AND BASIS POLYNOMIALS 


RETURN 


FIGURE 7 SUBROUTINE FINITE LOGIC FLOW 
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Table 3_ Central Processor Execution Time 
Comparison: Model 310 vs. 
Chebyshev Polynomial Model vs. 
Orthonormal Polynomial Model 
(average time per point). 


Chebyshev Orthonormal 


2.2 msec 2.5 msec 


dinates. The execution time ratios show that even for the worst cases 
the Orthonormal model and the Chebyshev model will be at least 40 times 
faster than an explicit point mass disturbance acceleration model such 
as Model 310. Note that by using the recurrence relations (24) we can 
reduce the execution times for Orthonormal so that they are almost 
identical to the ones produced by Chebyshev. 

A certain price must be paid to generate finite element fields. 
Execution times are shown in Table 4. The price per cell must be 
multiplied by the total number of cells in the finite element field and 


added to the price to generate the numbers f 6) Note that over 90% 


ij 
of the price paid per cell is just the generation of the gravity obser- 
vations from Model 310. Note the significant time-savings in using the 
orthonormal model. 

The total number of coefficients needed for various finite element 
field models of the 35° x 40° region covered by Model 310 are shown in 


Table 5. Table 5 assumes that Hoe will be constant at 300 km and that 


11 
300 km will be the maximum altitude to be modeled. For a given cell 
size, it shows how many finite element cells would be needed to model 
the region and how many total coefficients would be needed for various 
NORDERS to model the region. Based upon error analyses conducted on a 
small volume near the center of the 35° x 40° region, expected maximum 
error bounds were determined for certain cell sizes and NORDERS. The 
upper solid and dashed underline indicates an expected maximum absolute 
error in mgals. Solid underlines in both cases indicate that the ex- 
pected maximum absolute error was verified for small test volumes while 
dashes indicate the apparent trend. 


In order to study the maximum absolute error as functions of NORDER 


and cell size, it was decided to consider the maximum absolute error as 
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Table 4 


Central Processor Times to Create Finite Element Fields 


Time to generate matrix A 
(sec) 


Time to generate the 
observations using Model 
310 (per cell) 


Time to process matrix A 


Time to calculate the 
coefficients (per cell) 


Total Prices 


Time to generate and 
process matrix A 


Time to generate obser- 
vations and coefficients 
using Model 310 (per cell) 


32 


14.87 


1.24 


-43 


1.56 


15.30 


Chebyshev 


4 


-61 


29.27 


5.45 


1.43 


6.06 


30.70 


5 


23.55 


54.80 


80.29 


89.35 


41 


14.48 


58 


«22 


-99 


14.70 


Orthonormal 
4 5 6 


-70 1.57 3.16 


28.63 50.43 79.65 


-80 1.29 2.43 


-75 2.00 5.07 


1.50 2.86 5.59 


29.38 52.43 84.72 


Total number of Coefficients vs. N. vs. Cell Size 
For a 300 km x 35° x 40° Finite Element Field. 
(Chebyshev and Orthonormal) 


Cell Size 
Qoe11 . %e11) 


No. of Cells 


35x40 = 1400 ee 147,000 235,200 352,800 
\ 

28x32 = 896 53, 760 * 94,080 150,528 225,792 
\ \ 

24x27 = 648 38,880 68,040 \ 108,864 \ 163,296, 


\ 
20x23 = 460 27,600 48,300 77,280 \ 115,920 


\ 
18x20 = 360 21,600 37,800 60,480 \ 90,720 


—"——+ 
16x18 = 288 17, 280 30,240 48,384 72,576 
14x16 = 224 13,440 23,520 37,632 56,448 
13x15 = 195 11,700 20,475 32,760 49,140 


12x14 = 168 10,080 17,640 28,224 42,336 


Given: Hoel1 = constant = 300 km 


* = Finite element fields above this line are expected 
to have lerror| x < 1.5 mgals. 


+ = Finite element fields above this line are expected 
to have error| oy < 3.0 mgals. 


a function of cell size only for fixed NORDER (see Figure 8) and to 
consider maximum absolute error as a function of NORDER for fixed cell 
sizes (see Figure 9). For all cases, cell altitude oop was fixed at 
300 km. 

Figure 8 shows that the maximum absolute error varies approximately 
as the square of the cell dimension in A or $. For constant cell altitude 


and for small dee and $, the volume will vary approximately as the 


11 


square of he of > Hence the maximum absolute 


11 cell’ se reell % % cell" 
error is really proportional to cell volume. Figure 9 shows that the 
maximum absolute error varies approximately inversely as the square of 
NORDER. As the number of coefficients per cell varies approximately as 


the square of NORDER, the maximum absolute error may be considered to 


vary nearly inversely with the number of coefficients per cell. 


eee Genel 
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1x1 2x2 3x3 1x1 2x2 3x3 
Cell Size (degrees) Cell Size (degrees) 


c. NORDER = 5 d. NORDER = 6 


1x1 2x2 3x3 1x1 2x2 3x3 
Cell Size (degrees) Cell Size (degrees) 


Figure 8 Maximum Absolute Error versus Cell Size 
for Fixed NORDER's. 
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Figure 9 Maximum Absolute Error versus NORDER 
for Fixed Cell Sizes. 
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9.0 Example Finite Element Gravity Models 


Based upon the tradeoff studies conducted as described above, two 


example finite element gravity fields were generated. The 300 km x 10° x 
10° region with altitude 0 to 300 km, longitude 70°E to 80°E, and latitude 
25°S to 35°S was modeled. A finite element field with NORDER = 3 and 
with one hundred 300 km x 1° x 1° cells was generated as was a finite 
element field with NORDER = 5 and the same size and number of cells. An 
error analysis was conducted on both fields. The finite element model 

was compared on a point by point basis with Model 310 at sample points 
“away from the fit points" for altitudes of 1 m, 150 km, and 299 km. 

The error mean, the RMS error, and the maximum absolute error were 
determined for the three sets of sample points for both finite element 


fields and are shown in Tables 6 and 7. 


Table 6 


N=3 Example Finite Element Gravity Disturbance Field 
(Chebyshev and Orthonormal) 


H=lao 


ERROR 
RMS ERROR 


|ERROR| 44, 


H = 150 km 


ERROR 
RMS ERROR 


| ERROR| 4, 


H = 299 km 


ERROR 


RMS ERROR 


| ERROR| yay 


Table 7 


N=5 Example Finite Element Gravity Disturbance Field 
(Chebyshev and Orthonormal) 


H=l1n 


ERROR 
RMS ERROR 


|ERROR| 44 


H_= 150 km 
ERROR 
RMS ERROR 


| ERROR | yy 


H_ = 299 km 


ERROR 


RMS ERROR 


| ERROR| ax 
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10.0 Orthogonal Approach 


The main purpose of this section is to indicate a way to produce an 
accurate model for the geopotential itself. The Chebyshev and Ortho- 
normal approaches model the gravity disturbance acceleration compcnents 
independently and do not take into account the fundamental functional 


relationship 
g= W (32) 


Therefore, no rigorous representation for U can be obtained by using 
gravity components derived either by the Chebyshev or Orthonormal 
approach. 

A finite element procedure for modeling the geopotential that takes 
the functional relationships (32) between the geopotential and its 
gradient into account appears desirable. Not only will this enable us 
to obtain an approximation for U, but a significant reduction in the 
overall cost of determining gradient models is anticipated. 

Let us consider the following preliminary approximation of the geo- 


potential disturbance function 


i ip 

6U = c f » k= n-i-j (33) 

n=0 i=0 j=0 13K idk 

where “45k are modeling constants to be determined, and £55 (ee¥ 92) are 


arbitrary basis functions. Using Eq. (32) and (33) we can write the 


acceleration disturbances as follows 


ies : 
6g, = c f (34a) 
x n=0 4=0 j=0 ijk ijk 
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seeaennaremammaapaamssasasassssmsmsssssmsmssssssmrpmmmmsssill 


——S ae <a eT 


N on i-n y 


ég. = c f (34b) 
87” ko 120 geo RY 
La : 
6g. = c f (34c) 
* n=0 izo j=o TK 45K 
with 
* 9 y 3 " 3 
Sage 7” Ge f2yu* “eqn ~ By 19k? Taye * Se £198 (35) 


The local modeling coefficients c are to be determined via a least 


ijk 
Squares approximation technique, taking into account the local observa- 
tions of all three gravity components (68,68, .68,) simultaneously. 


The normal equations can be written in the usual manner 


c= (Atay asp” (36) 
where 
“0007 S8%(1) 
“001 6g*(1) 
rs y 
c= c 6 = (37) 
£ o10. f-* “& 5g%(1) 
“N00 Sex (m) 
bpe% 
By (m) 
5% (m) 


are the column matrix of modeling constants c and the column matrix 


ijk 


of gravity disturbance observations, respectively. Furthermore, 


————EE 


(m) aie 


ooo™  — F9914™ ~~ Fg 


y y y Mf 
foo08™ fog = Eggg(™ =— ees Egg (™ 


Zz z z z 
|foo0™ £991 (™ £o19™ eae £00 ™ 


in which m indicates the total number of local measurements of the dis- 
turbance acceleration. 


It is now possible to evaluate the constants c from Eq. (36) in 


ijk 


the same manner as the Chebyshev approach. However, we wish to introduce 


a significant simplification by using orthogonal basis functions, similar 


to the Orthonormal approach. To this end, let us write the expression 


for a generic element of the matrix Ala as follows 


m 
dF (2) (39) 
£=1 ijk,apy 
where 
x x y y z z 
Pade = fii) fay 6) + fii Egy 6%) + £556) £ py \*) (40) 


Let us now introduce basis functions of the form 


£5 5 (%o¥ 92) = £, (x) £,(y) f(z) (41) 
so that, 
y Zz 
tay = fif Sf, > Fisk = a ot ; fii = hae be (42) 
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where the derivatives denoted by primes are taken with respect to the 


argument. It is clear that the generic element (39) can be written as 
, 1 
yet. IE, £, Ye.£,, 
x y z 
' 8 
+ if, If, f, yet, (43) 
x y z 
tet 
+ Tee If ,f, pete’ 
x y z 


where the respective summations run over all the observations in the 
respective directions. 


Next, let us assume that the sets of functions 


EY, ey th) 


j 


are orthonormal over the interval of validity and that similarly 


HE iy (Ew AEE 


are orthogonal with respect to the same interval, then it is clear that 


ah = D becomes diagonal with diagonal elements Pi dk given by 


iz 12 12 2 
Disk BES + fs + yey 
x y z 
Equation (36) becomes 


“ * 
c#D 1,7 6g 
or, more explicitly, 


Rae. Sha a ee Sot yt aE (say 
“age” "ag. gh gk Me,” Teak, 8 Tage, FGe)e ) n 
The only remaining task is to find an orthonormal set of functions 


{£, (x) } such that the derivative set {£, (x) } is also orthogonal over the 


same interval of validity. Note that the same set can be used in all 
three directions x,y,z. 

To date, only NORDER = 3 has been tried. As expected, however, not 
enough accuracy is obtained and further work is underway for higher 
NORDERs. It is expected that the number of modeling constants Cs 5k 
will be significantly reduced compared to the total number of modeling 
constants ae eee Crk? necessary to model the gravity disturbance 
components independently. This fact, not only increases the computa- 


tional speed of the algorithm, but also reduces the storage requirements. 


Moreover, it is now possible to obtain a valid approximation for the 


geopotential itself. 
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12.0 Software Listings 


Page 
Point Mass Model 310 58 
| Program MASPOS 59 
Listing of Model 310 Point Mass Coordinates and GM Values 62 
Chebyshev Polynomial Finite Element Software 83 
Program LOCALG 84 
i Subroutine CHEBY 91 
Subroutine PTMASS 92 
Program FINEG 95 
Subroutine ALSQ 98 
Program FINTES 102 
Subroutine FINITE 108 
Listing of Chebyshev Finite Element Disturbance Gravity 
Acceleration Output 116 
Orthornormal Polynomial Finite Element Software 117 
Program LOCALG* ; 118 
Subroutine ORTHO 125 
Subroutine MULT 127 
Subroutine PTMASS* 128 
Program ORTHOG 131 
Subroutine ALSQI* 134 
Program FINTES* 136 
Subroutine FINITE* 142 
Listing of Orthogonal Finite Element Disturbance Gravity 
° Acceleration Output** 149 


x 
These Orthonormal programs and subroutines represent minor modifi- 
cations of the corresponding Chebyshev Polynomial Software. 


ak 
The numerical results (for fixed NORDERS) are identical for the two | 


approaches. The orthonormal vs. Chebyshev tradeoffs are with regard to 
efficiency, not accuracy, see the above discussion in Section 8.0. 
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